function [pf,nrm] = disclegfunc(f,n)
% this function takes in a discrete function f
% and computes the nth order discrete legendre 
% approximation to f and the corresponding L2 norm

N = length(f)-1;      % max degree
k = 0:N;            % interval
P = zeros(N,N+1);
P(1,1:N+1) = 1; 
P(2,1:N+1) = 1 - (2*k)/N;

for i = 1:n-1 %otherwise N-1
    P(i+2,1:N+1) = ( (2*i+1)*(N-2*k).*P(i+1,:) - i*(N+i+1)*P(i,:) ) ./ ...
        ((i+1)*(N-i)) ;
end

for i = 0:n
    a = (2*i+1);
    b = factorial(N+i+1)/factorial(N);
    c = factorial(N)/factorial(N-i);
    w(i+1)=a*c/b;
end
ff=zeros(1,N+1);
for ord=1:n+1
    ff(ord) = w(ord) * f * P(ord,:)';
end

pf = zeros(1,N+1);
for ord=1:n+1
    pf = pf + ff(ord)*P(ord,:);
end
ff'
plot( k, f, 'k', k, pf, 'r--')

nrm = norm(f-pf);
